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We present new techniqes for evolving binary black hole systems which allow the accurate determi- 
nation of gravitational waveforms directly from the wave zone region of the numerical simulations. 

Rather than excising the black hole interiors, our approach follows the “puncture” treatment of 
black holes, but utilzing a new gauge condition which allows the black holes to move successfully 
through the computational domain. We apply these techniques to an inspiraling binary, model- 
ing the radiation generated during the final plunge and ringdown. We demonstrate convergence of 
the waveforms and and good conservation of mass-energy, with just over 3% of the system’s mass 
converted to gravitional radiation. 
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Coalescing comparable mass black hole binaries are 
prodigious sources of gravitational waves. The final 
merger of these systems, in which the black holes leave 
their quasicircular orbits and plunge together to produce 
a highly distorted black hole that “rings down” to a qui- 
escent Kerr state, will produce a strong burst of gravita- 
tional radiation. Such mergers are expected to be among 
the strongest sources for ground-based gravitational wave 
detectors, which will observe the mergers of stellar-mass 
and intermediate mass black hole binaries, and the space- 
based LISA, which will detect mergers of massive black 
hole binaries. Observations of these systems will provide 
an unprecedented look into the strong-field dynamical 
regime of general relativity. With the first-generation 
of ground-based interferometers reaching maturity and 
LISA moving forward through the formulation phase, the 
need for accurate merger waveforms has become urgent. 

Such waveforms can only be obtained through 3-D nu- 
merical relativity simulations of the full Einstein equa- 
tions. While this has proven to be a very challenging 
undertaking, new developments allow an optimistic out- 
look. Full 3-D evolutions of binary black holes, in which 
regions within the horizons have been excised from the 
computational grid, have recently been carried out. Us- 
ing co-rotating coordinates, so that the holes remain fixed 
on the grid as the system evolves, a binary has been 
evolved through a little more than a full orbit [1] as well 
as through a plunge, merger, and ringdown [2], though 
without being able to extract gravitational waveforms. 
More recently, a simulation in which excised black holes 
move through the grid in a single plunge-orbit, merger, 
and ringdown has been accomplished, with the calcula- 
tion of a waveform [3]. 

In this Letter , we report the results of new simula- 
tions of inspiraling binary black holes through merger 
and ringdown. These have been carried out using new 
techniques which allow the black holes to move through 
the coordinate grid without the need for excision [17]. 
Using fixed mesh refinement, we are able to resolve both 
the dynamical region where the black holes inspiral (with 


length scales ~ M, where M is the total system mass) 
and the outer regions where the gravitational waves prop- 
agate (length scales ~ (10 - 100)M). Using an outer 
boundary at 128 AT, we evolve the system to well beyond 
t « 100 M, extract gravitational waveforms and demon- 
strate that they are 2 nd -order convergent. 

We start by setting up “puncture” initial data for 
equal mass binary black holes [4]. The metric on the 
initial spacelike slice takes the form — -0 4 Sij , where 
i,j = 1,2, 3, and the conformal factor 0 = 0 bl + u. The 
static, singular part of the conformal factor has the form 
0 bl = 1 + ELi m ni 2|r — f n |, where the n th black hole 
has mass m n and its located at r n . The nonsingular func- 
tion u is obtained by solving the Hamiltonian constraint 
equation using AMRMG [5] . We use parameters so that the 
black holes have proper separation 4.99 M, and the sys- 
tem has total mass M — 1.008 and angular momentum 
J = 0.779 M 2 . This corresponds to the run QC0 studied 
in Ref. [6]. 

In the standard puncture implementation, 0 bl Is fac- 
tored out and handled analytically; only the regular parts 
of the metric are evolved. In this case, the punctures 
remain fixed on the grid while the binary evolves. How- 
ever, the stretching of the coordinate system that ensues 
is problematical. First, as the physical distance between 
the black holes shrinks, certain components of the met- 
ric must approach zero, causing other quantities to grow 
uncontrollably. Second, a corotating coordinate frame 
(implemented by an appropriate angular shift vector) is 
necessary to keep the orbiting punctures fixed on the grid; 
this causes extremely superluminal coordinate speeds at 
large distances from the black holes and, in the case of a 
Cartesian grid, incoming noise from the outer boundary. 

Our approach is to allow the punctures to move freely 
through the grid, by not factoring out the singular part 
of the conformal factor but rather evolving it insepara- 
bly from the regular part. Initially, we follow the stan- 
dard puncture technique and set up the binary so that 
the centers of the black holes are not located at a grid 
point. Taking numerical derivatives of -0 bl effectively 
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regularizes the puncture singularity using the smoothing 
inherent in the finite differences. These regularized data 
are then evolved numerically. Since the centers of the 
black holes remain in the z ~ 0 plane, they do not pass 
through gridpoints in our cell-centered implementation. 

We evolve this data with the Hahndol code, which uses 
a conformal formulation of Einstein’s evolution equations 
on a cell-centered numerical grid [7] with a box-in-box 
resolution structure implemented via Paramesh [8] The 
innermost refinement region is a cube stretching from 
—2 M to 2 M in all 3 dimensions, and has the finest res- 
olution hf. The punctures are placed within this region 
on the y — axis in the z = 0 plane; we impose equatorial 
symmetry throughout. We performed three simulations 
with identical grid structures, but with uniformly differ- 
ing resolutions. In the most refined cubical region the 
resolutions were hf = M/16, M/24, and M/32. Subse- 
quent boxes of doubled size have half the resolution. We 
use 8 boxes to put the outer boundary at 128M, causally 
disconnected from the wave extraction region through 
most of the run. We use 4th order finite differencing 
for the spatial derivatives except for the advection of the 
shift, which is performed with 2 nd -order, mesh-adapted 
differencing [9], and we use 2 nd -order time stepping via 
a three-step iterative Crank-Nicholson scheme. 

In our new approach, the free evolution of punctures is 
made possible by a modified version of a common coor- 
dinate condition known as the Gamma- freezing shift vec- 
tor, which drives the coordinates towards quiescence as 
the merged remnant black hole also becomes physically 
quiescent. Our modified version retains this “freezing” 
property, yet is suitable for motile punctures. Specifically 
we use dt/3 1 = \aB l and d t B x = d t f 1 — ft djT 1 — 
which incorporates two critical changes to the standard 
Gamma- freezing condition. A factor -0 bl of the confor- 
mal factor, originally used to ensure that the shift van- 
ishes at the puncture, has been removed in order to allow 
the punctures to move. Also, a new term has been added 
(— ft djY % ) which facilitates more stable and accurate evo- 
lution of moving punctures by eliminating a zero-speed 
mode (which was otherwise found to create a “punc- 
ture memory” effect as errors grew in place [10]). Along 
with this shift condition, we use the standard singularity- 
avoiding, 14-log slicing condition on the lapse. 

Fig. 1 shows the error in the Hamiltonian constraint 
Ch at 2 different times. The peak violation near the 
puntures does not leak-out, or grow with time, but stays 
well- confined even though the punctures and horizons 
move across the grid. Overall, we get 2 nd — order conver- 
gence away from the horizons to well beyond the wave 
extraction region for the entire course of the run. There 
is no indication of exponentially growing constraint vio- 
lations which have plagued many numerical simulations 
with black holes fixed in place on the numerical grid. 

One way to get a picture of the motion of the black 
holes is to look at the location of the black hole appar- 



FIG. 1: Hamiltonian constraint error Ch for hf *= M/24 
and M/32, at two times when a puncture is near to crossing 
the positive x-axis. The data are scaled such that the lines 
should superpose in the case of perfect 2 nd -order convergence. 
The inset shows that Ch is well-behaved in the region near 
the punctures. The horizontal lines indicate the approximate 
location of the apparent horizons; at the later time a common 
horizon has formed. 



FIG. 2: The positions of the apparent horizons at times t = 
0,5, 10, 15, and 20M for our M/16 run. The curve shows the 
trajectories of centroids of the individual apparent horizons. 

ent horizons at different times. Fig. 2 shows the loca- 
tions of a sequence of apparent horizons (calculated us- 
ing the AHFINDERDIRECT code[ll]) where they cross the 
x-y- plane for our hf = M/16 run. In the coordinates of 
our simulation, the black holes undergo about one-half 
orbit before forming a common horizon. 

We extract the gravitational waves generated by the 
merger using the technique explained in detail in [12]. 
Fig. 3 shows the dominant l = 2,m = 2 components 
of the Weyl curvature scalar extracted at 2 different 
radii from the medium and high resolution runs. For 
each resolution, the time-shifted and rescaled waveforms 
computed at different extraction radii are nearly indistin- 
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FIG. 3: Real part of r4>4 extracted from the numerical simu- 
lation on spheres of radii tex — 20, and 40 M for the medium 
and high resolution runs. The waveforms extracted at differ- 
ent radii have been rescaled by l /vex and shifted in time to 
account for the wave propagation time between the extrac- 
tion spheres. At high resolution (hf = M/32) there is no 
discernible dependence on extraction radius. For comparison, 
we show Lazarus waveforms from Ref. [6]. 

guisable, indicating that the waves travel cleanly across 
refinement boundaries and have the expected 1/r falloff. 
In addition, the two highest resolution waveforms differ 
only by a slight phase shift, and a by few percent in am- 
plitude. For comparison we have also included the QCO 
Lazarus waveforms from Sec. V of Ref. [6] . These were 
extracted by approximately matching the later portion of 
brief numerical simulations onto a perturbed black hole 
[13, 14] at transition time 10 M. 

Fig. 4 shows the convergence of the extracted waves 
throughout the run. The difference between the two high- 
est resolutions is roughly 90 degrees out of phase with 
the waveform, corresonding to a small phase-shift in the 
waveform, possibly caused by a small difference in the 
orbital trajectories. 

The gravitational waveforms also contain physical in- 
formation about the radiation, including the energy E 
and angular momentum J carried away by the radiation. 
We calculate dE/dt and dJ/dt from time integrals of all 
l = 2 waveform components using Eqs. (5.1) and (5.2) in 
[6]. Integrating dE/dt gives the energy loss as a function 
of time; this should be the same as M — M A dm, where 
Madm is the ADM mass extracted at on a sphere of 
sufficiently large radius [15]. Fig. 5 shows a comparison 
between M~E and an independent calculation of Madm 
at two extraction radii tex = 40 and 50 M. The striking 
consistency between the two calculations as the radiation 
passes indicates good energy conservation in the simula- 
tion. Shortly before the arrival of the radiation the ADM 
mass measurement is affected by a transient non-physical 
pulse in the gauge evolution, though the pulse does not 
affect the radiation measurement. 
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FIG. 4: Differences of the real part of r^ 4 for resolutions of 
hf = M/16, M/24, and M/32 appropriately scaled such that 
for perfect 2 nd -order convergence the lines would lay on top 
of each other. 



FIG. 5: Conservation of mass-energy for the highest resolu- 
tion case, hf — M/32 . We compare the ADM mass Madm 
with the mass remaining, M — E, after gravitational radia- 
tion energy loss E . The good agreement, based on extraction 
spheres at vex = 40 and 50 M, indicates conservation of en- 
ergy in the simulation. 


The total radiated energy calculated from the wave- 
forms extracted at Tex = 20, 30, 40 and 50M in the high- 
est resolution run has the values E/M —0.0304, 0.0312, 
0.0317 and 0.0319, respectively. While these values vary 
significantly with tex (even extracting at these relatively 
large radii), they are neatly consistent with a 1/tex 
falloff to an asymptotic value of 0.0330 with an uncer- 
tainty in the extrapolation of < 1%. In Table I we give 
the total radiated energies and angular momenta extrap- 
olated as tex oo. For comparison we also include the 
Lazarus values, as well as values from the AEI group [2] 
which did not determine waveforms, but estimated the 
radiative losses based on the state of the final black hole 
horizon in runs including the QCO case. Our lowest res- 
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M/16 

M/ 24 

M/32 

Lazarus AEI 

E/M 

0.0516 

0.0342 

0.0330 

0.025 

0.030 

J/M 2 

0.208 

0.140 

0.138 

0.10 

0.17 


TABLE I: The radiated energy E and angular momentum J 
carried away by gravitational radiation in our simulations. 
Our values are comparable with earlier estimates from the 
A El (via horizon analysis) and Lazarus via perturbation tech- 
niques. 
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olution run clearly over-estimates the radiation energy 
and angular momentum while our higher resolution re- 
sults are in closer agreement with the AEI value for the 
energy, and close to the 20% level of confidence suggested 
in Ref. [6]. 

In conclusion, we have calculated a gravitational ra- 
dition waveform directly via numerical simulation of an 
inspiraling configuration of binary black holes. Starting 
with equal mass puncture black holes and parameters 
corresponding to the QCO model of Ref. [6], we release 
the holes and allow them to move through the grid with- 
out excision. A new gauge condition is used that allows 
us to accurately evolve this system from the initial inspi- 
ral orbit through merger and ringdown. The simulations 
converge with increasing resolution to 2 nd -order, leading 
to a 2 nd -order convergence of the waveform. These wave- 
forms have the correct 1 jr falloff and agree to a great 
extent with approximatively calculated ones. Our sim- 
ulations show good energy conservation as indicated by 
comparing the change in ADM mass with the radiated 
energy. The QCO configuration provides a model for the 
final plunge of the two black holes and the subsequent 
ringdown. In this brief burst of gravitational radiation 
we find that just over 3% of the system’s initial mass- 
energy is carried away in gravitational waves. 

The new gauge allows simulations to remain accurate 
far longer than previous standard puncture techniques. 
Our treatment will generalize, allowing us to study ra- 
diation generation in simulations of a variety of initial 
black hole configurations. Using adaptive mesh refine- 
ment, we plan to apply these techniques to study binaries 
beginning from larger initial separation, which are ex- 
pected to provide more realistic models corresponding to 
astrophysical systems. For further understanding of such 
model dependence, we will compare results from simula- 
tions beginning with different initial data models. We 
will also study the effects of unequal black hole masses, 
and the individual black hole spins. 

We thank David Brown for providing AMRMG; and 
Jonathan Thornburg for providing AHF I NDERD I RECT , 
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